
#include "VortexTransportProblem.h"

template<>
InputParameters validParams<VortexTransportProblem>()
{
  InputParameters params = validParams<EulerProblem>();
  params.addParam<Real>("R_gas",  287.15, "气体常数");
  params.addRequiredParam<Real>("strength", "气体常数");
  return params;
}

VortexTransportProblem::VortexTransportProblem(const InputParameters &params) :
	EulerProblem(params),
	_R_gas(getParam<Real>("R_gas")),
	_cp(_gamma/(_gamma-1)*_R_gas),
	_center(0.05, 0.05, 0),
	_radius(0.005),
	_strength(getParam<Real>("strength")),
	_pressure_inf(1E+05),
	_temperature_inf(300),
	_density_inf(_pressure_inf/_R_gas/_temperature_inf),
	_velocity_inf(_mach*sqrt(_gamma*_R_gas*_temperature_inf))
{

}


void VortexTransportProblem::valueExact(Real* value, Real t, const Point& p)
{
	Real x = p(0);
	Real y = p(1);
	Real z = p(2);
	Real xc = _center(0);
	Real yc = _center(1);
	Real r2 = (p-_center).norm_sq()/_radius/_radius;

	Real du = -(y-yc)*exp(-r2/2.)*(_velocity_inf*_strength/_radius);
	Real dv =  (y-yc)*exp(-r2/2.)*(_velocity_inf*_strength/_radius);
	Real dT = exp(-r2)*(_velocity_inf*_strength)*(_velocity_inf*_strength)/_cp/2.;

	Real u = _velocity_inf*(1+du), v = _velocity_inf*dv;
	Real T = _temperature_inf - dT;
	if(fabs(dT) > 0.0001) std::cout << dT << std::endl;
	Real density = _density_inf*pow(T/_temperature_inf, 1/(_gamma-1));
	Real pre = density*_R_gas*T;
	RealVectorValue mom(density*u, density*v, 0);

	value[0] = density;
	value[1] = mom(0);
	value[2] = mom(1);
	value[3] = mom(2);
	value[4] = pre/(_gamma-1) + 0.5*mom.norm_sq()/density;
}

void VortexTransportProblem::initialCondition(Real* value, Real t, const Point& p)
{
	valueExact(value, t, p);
}
